import numpy as np
from matplotlib import image
import matplotlib.pyplot as plt
import matplotlib.colors as clr
from scipy.fftpack import dct
from scipy.fftpack import idct
import pandas as pd
import PIL as p
from matplotlib.colors import Normalize
%matplotlib inline
"""
Returns an array of type:
- [x, y] for grayscale images
- [x, y, [R, G, B]] for RGB images
- [x, y, [R, G, B, A]] for RGBA images
"""
def read_image(filename):
img = image.imread(filename)
return img
"""
Funcao para criar um colormap;
Parametros:
1.cmin-> tuplo com valores (r,g,b) minimos
2.cmax-> tuplo com valores (r,g,b) maximos
Devolve o objeto correspondente ao colormap
"""
def create_colormap(cmin : tuple[float], cmax: tuple[float]):
return clr.LinearSegmentedColormap.from_list('', [cmin, cmax], N=256)
"""
Função para mostrar uma imagem. Aceita um colormap definido pelo utilizador ou os do matplotlib
"""
def show_image(img, colormap = None):
plt.figure(figsize=(8,8))
# Imagens com apenas uma coponenete: R, G, B ou Grayscale
if len(img.shape) == 2:
plt.imshow(img, cmap = colormap)
else:
if colormap != None:
new_img = img[:, :, 0]
plt.imshow(new_img, cmap = colormap)
else:
plt.imshow(img)
plt.axis('off')
plt.show()
#img = read_image('imagens/peppers.bmp')
img = read_image('imagens/barn_mountains.bmp')
#img = read_image('imagens/logo.bmp')
show_image(img)
| Barn | Peppers | Logo | |
|---|---|---|---|
| size | 356.5KB | 589.9KB | 421.6KB |
Foi utilizado o Adobe Photoshop 2019 para comprimir as imagens.
| Barn | Peppers | Logo | |
|---|---|---|---|
| Low | 43.4KB | 35.2KB | 21.9KB |
| Medium | 51.5KB | 41.3KB | 23.1KB |
| High | 67.5KB | 57.7KB | 27.3KB |
| Barn | Peppers | Logo | |
|---|---|---|---|
| Low | 8.2:1 | 16.7:1 | 19.3:1 |
| Medium | 6.9:1 | 14.3:1 | 18.3:1 |
| High | 5.2:1 | 10.2:1 | 15.4:1 |
| Barn | Peppers | Logo | |
|---|---|---|---|
| Low | não muito evidente | muito evidente | evidente |
| Medium | nada evidente | não muito evidente | evidente |
| High | nada evidente | nada evidente | nada evidente |
A compressão é maior em imagens cujos pixels adjacentes têm uma cor semelhante.
Isto sugere que o codec do Jpeg utiliza compressão baseada em adjacência.
Visualmente, notam-se maiores diferenças em locais onde há alterações bruscas na cor do pixel.
low = image.imread('imagens/Low/peppers.jpg')
low.shape
(384, 512, 3)
medium = image.imread('imagens/Medium/peppers.jpg')
medium.shape
(384, 512, 3)
# ColorMaps
cm_gray = clr.LinearSegmentedColormap.from_list('gray', [(0,0,0), (1, 1, 1)], N = 256)
cm_red = clr.LinearSegmentedColormap.from_list('red', [(0,0,0), (1, 0, 0)], N = 256)
cm_green = clr.LinearSegmentedColormap.from_list('green', [(0,0,0), (0, 1, 0)], N = 256)
cm_blue = clr.LinearSegmentedColormap.from_list('blue', [(0,0,0), (0, 0, 1)], N = 256)
# Adaptado de https://jakevdp.github.io/PythonDataScienceHandbook/04.07-customizing-colorbars.html
cmap = plt.cm.get_cmap(cm_gray)
gray = cmap(np.arange(cmap.N))
cmap = plt.cm.get_cmap(cm_red)
red = cmap(np.arange(cmap.N))
cmap = plt.cm.get_cmap(cm_green)
green = cmap(np.arange(cmap.N))
cmap = plt.cm.get_cmap(cm_blue)
blue = cmap(np.arange(cmap.N))
fig, ax = plt.subplots(4, figsize=(10, 5),subplot_kw=dict(xticks=[], yticks=[]))
ax[0].imshow([gray], extent=[0, 10, 0, 1])
ax[1].imshow([red], extent=[0, 10, 0, 1])
ax[2].imshow([green], extent=[0, 10, 0, 1])
ax[3].imshow([blue], extent=[0, 10, 0, 1])
<matplotlib.image.AxesImage at 0x7ff7732de700>
"""
Separar uma imagem RGB nos seus componentes
"""
def separate_components(img):
r = img[:, :, 0]
g = img[:, :, 1]
b = img[:, :, 2]
return r, g, b
r, g, b = separate_components(img)
show_image(r, cm_red)
show_image(g, cm_green)
show_image(b, cm_blue)
"""
Juntar as coponentes R, G e B para formar uma imagem
"""
def join_components(r, g, b):
return np.dstack((r, g, b))
img = join_components(r, g, b)
show_image(img)
"""
Recebe uma imagem e altera as suas dimensões (m,n) para (16*p, 16*q).
Isto é realizado através da cópia da ultima coluna/linha até atingir o valor multiplo de 16.
Conta o numero de linhas/colunas adicionadas, (x,y) : 0<=x,y<=15.
Devolve um tuplo com as dimensoes incrementadas (x,y) e a imagem com as novas dimensoes (16*p, 16*q)
"""
def padding(img : np.array):
img = img.copy()
shape = img.shape
x,y = (16-img.shape[0]%16)%16, (16-img.shape[1]%16)%16
h_padding = np.repeat(img[-1:,:,:], x, axis = 0)
img = np.concatenate((img, h_padding), axis = 0)
v_padding = np.repeat(img[:,-1:,:], y, axis = 1)
img = np.concatenate((img, v_padding), axis = 1)
return shape, img
"""
Recebe uma imagem e as dimensões originais dela.
Faz slice da imagem, elimando todos os elementos incrementados no padding;
Devolve a imagem original
"""
def unpad(img : np.array, shape):
img = img.copy()
x,y,z = shape
img = img[:x, :y]
return img
show_image(img)
shape, pad_img = padding(img)
show_image(pad_img)
unpad_img = unpad(img, shape)
show_image(unpad_img)
# Verificação que fica igual ao original
print(np.all(np.equal(img, unpad_img)))
print(f"Original shape: {img.shape}")
print(f"Padding shape: {pad_img.shape}")
print(f"After Padding and Unpadding shape: {unpad_img.shape}")
True Original shape: (297, 400, 3) Padding shape: (304, 400, 3) After Padding and Unpadding shape: (297, 400, 3)
"""
Converte imagem no formato RGB para imagem no formato yCbCr;
"""
def rgb_to_ycbcr(img : np.array):
img.copy().astype(float)
y_cb_cr_mat = np.array([ [0.299 , 0.587 , 0.114 ]
,[-0.168736, -0.331264, 0.5 ]
,[0.5 , -0.418688, -0.081312] ], dtype=float)
y = y_cb_cr_mat[0,0] * img[:,:,0] + y_cb_cr_mat[0,1] * img[:,:,1] + y_cb_cr_mat[0,2]*img[:,:,2]
cb = y_cb_cr_mat[1,0] * img[:,:,0] + y_cb_cr_mat[1,1] * img[:,:,1] + y_cb_cr_mat[1,2]*img[:,:,2] + 128
cr = y_cb_cr_mat[2,0] * img[:,:,0] + y_cb_cr_mat[2,1] * img[:,:,1] + y_cb_cr_mat[2,2]*img[:,:,2] + 128
y_cb_cr = np.dstack((y, cb, cr))
return y_cb_cr
"""
Converte imagem no formato YCbCr para imagem no formato RGB
"""
def ycbcr_to_rgb(img : np.array):
img = img.copy().astype(float)
y_cb_cr_mat_inv = np.linalg.inv(
np.array([ [0.299 , 0.587 , 0.114 ]
, [-0.168736, -0.331264, 0.5 ]
, [0.5 , -0.418688, -0.081312] ], dtype=float)
)
y = img[:, :, 0]
cb = img[:, :, 1] - 128
cr = img[:, :, 2] - 128
r = y + y_cb_cr_mat_inv[0,2]*cr
g = y + y_cb_cr_mat_inv[1,1]*cb + y_cb_cr_mat_inv[1,2]*cr
b = y + y_cb_cr_mat_inv[2,1]*cb
rgb = np.dstack((r,g,b))
rgb = np.round(rgb)
rgb[rgb > 255] = 255
rgb[rgb < 0] = 0
return rgb.astype(np.uint8)
ycbcr = rgb_to_ycbcr(pad_img)
rgb = ycbcr_to_rgb(ycbcr)
# Verificação que fica igual ao original
#print(np.all(np.equal(img, rgb)))
show_image(img)
show_image(rgb)
# Método do colormap para representar os canais Cb e Cr:
# https://stackoverflow.com/questions/28638848/displaying-y-cb-and-cr-components-in-matlab
cm_cb = clr.LinearSegmentedColormap.from_list('cb', [(0.5, 0.5, 0), (0.5, 0.5, 1)], N = 256)
cm_cr = clr.LinearSegmentedColormap.from_list('cr', [(0, 0.5, 0.5), (1, 0.5, 0.5)], N = 256)
cmap = plt.cm.get_cmap(cm_cb)
cm_cb_rep = cmap(np.arange(cmap.N))
cmap = plt.cm.get_cmap(cm_cr)
cm_cr_rep = cmap(np.arange(cmap.N))
fig, ax = plt.subplots(2, figsize=(10, 2),subplot_kw=dict(xticks=[], yticks=[]))
ax[0].imshow([cm_cb_rep], extent=[0, 10, 0, 1])
ax[1].imshow([cm_cr_rep], extent=[0, 10, 0, 1])
<matplotlib.image.AxesImage at 0x7ff7718fa2b0>
y, cb, cr = separate_components(ycbcr)
show_image(y, cm_gray)
show_image(cb, cm_cb)
show_image(cr, cm_cr)
f, ax = plt.subplots(nrows=3,ncols=2,figsize=(12,12))
ax[0,0].imshow(y, cm_gray)
ax[0,1].imshow(r, cm_red)
ax[1,0].imshow(cb, cm_cb)
ax[1,1].imshow(g, cm_green)
ax[2,0].imshow(cr, cm_cr)
ax[2,1].imshow(b, cm_blue)
for i in ax.flatten():
i.set_axis_off()
Verifica-se que o canal Y da imagem tem muito mais definição do que o Cb e do que o Cr.
Em relação aos canais R,G e B, os canais R e G são os que apresentam maior semelhança visual com o canal Y. O canal G é o que tem mais peso no calculo do canal Y, seguido do R e finalmente do B.
Isto deve-se ao facto dos niveis de luminância ser maior precisamente nos canais G, seguido do R, seguido do B.
Com isto, concentra-se a maior parte da informação da luminância num dos canais, e informação sobre a crominância, que é menos importante para o olho humano nos outros (Cb e Cr).
Isto permite que haja maior compressão nos canais Cb e Cr, pois esta informação afeta menos a perceção do olho humano.
Semana 2¶
Aplicação dos primeiros passos destrutivos do CODEC jpeg
- Downsample com os rácios 4:2:2 e 4:2:0
- DCT aplicada á imagem toda
"""
(Y:Cb:Cr)
(4, 2, 0) -> São removidos pixeis nas linhas e colunas
(4, 2, 2) -> São removidos pixeis apenas nas colunas
"""
def sub_sample(y, cb, cr, downsample_ratio):
y = y.copy()
cb = cb.copy().astype(float)
cr = cr.copy().astype(float)
if downsample_ratio == (4,2,0):
ratio = round(downsample_ratio[0]/downsample_ratio[1])
cb = cb[::ratio,::ratio]
cr = cr[::ratio,::ratio]
else:
v_ratio = round(downsample_ratio[0]/downsample_ratio[1])
h_ratio = round(downsample_ratio[0]/downsample_ratio[2])
cb = cb[:, ::v_ratio]
cr = cr[:, ::h_ratio]
return y,cb,cr
"""
funciona em -- (4,2,0) e (4,2,2) --
Parametros:
1. y --------------> matriz com os valores da luminancia Y
2. cb -------------> matriz downsampled de cb
3. cr -------------> matriz downsample de cr
4. upsample_ratio -> tuplo com 3 inteiros correspondentes ao ratio de downsample (Y:Cb:Cr)
"""
def const_up_sample(y, cb, cr, upsample_ratio):
y = y.copy()
cb = cb.copy()
cr = cr.copy()
cbArr = np.zeros(shape = y.shape, dtype=float)
crArr = np.zeros(shape = y.shape, dtype=float)
if upsample_ratio[-1] == 0:
ratio = int(upsample_ratio[0]/upsample_ratio[1])
cbArr[::ratio,::ratio] = cb
cbArr[1::ratio,::ratio] = cb
cbArr[::ratio,1::ratio] = cb
cbArr[1::ratio,1::ratio] = cb
crArr[::ratio,::ratio] = cr
crArr[1::ratio,::ratio] = cr
crArr[::ratio, 1::ratio] = cr
crArr[1::ratio,1::ratio] = cr
else:
cb_ratio = int(upsample_ratio[0]/upsample_ratio[1])
cr_ratio = int(upsample_ratio[0]/upsample_ratio[2])
cbArr[:,::cb_ratio] = cb
cbArr[:,1::cb_ratio] = cb
crArr[:,::cr_ratio] = cr
crArr[:,1::cr_ratio] = cr
return y,cbArr,crArr
"""
funciona em -- (4,2,0) e (4,2,2) --
Parametros:
1. y --------------> matriz com os valores da luminancia Y
2. cb -------------> matriz downsampled de cb
3. cr -------------> matriz downsample de cr
4. upsample_ratio -> tuplo com 3 inteiros correspondentes ao ratio de downsample (Y:Cb:Cr)
"""
def linear_up_sample(y, cb, cr, upsample_ratio):
y = y.copy()
cb = cb.copy()
cr = cr.copy()
new_cb = np.zeros(shape = y.shape, dtype=float)
new_cr = np.zeros(shape = y.shape, dtype=float)
if upsample_ratio[2] == 0:
ratio = int(upsample_ratio[0]/upsample_ratio[1])
#--------------------Cb interpolation----
cb_cols_mean = (cb[:,:-1] + np.roll(cb, -1, 1)[:,:-1])//2
cb_cols_mean = np.concatenate((cb_cols_mean, cb[:,-1:]), axis = 1)
new_cb[::2,::2] = cb
new_cb[::2,1::2] = cb_cols_mean
cb = new_cb[::2]
cb_lines_mean = (cb[:-1] + np.roll(cb, -1, 0)[:-1])//2
cb_lines_mean = np.concatenate((cb_lines_mean, cb[-1:]), axis = 0)
new_cb[1::2,:] = cb_lines_mean
#------------------------Cr interpolation----
cr_cols_mean = (cr[:,:-1] + np.roll(cr, -1, 1)[:,:-1])//2
cr_cols_mean = np.concatenate((cr_cols_mean, cr[:,-1:]), axis = 1)
new_cr[::2,::2] = cr
new_cr[::2,1::2] = cr_cols_mean
cr = new_cr[::2]
cr_lines_mean = (cr[:-1] + np.roll(cr, -1, 0)[:-1])//2
cr_lines_mean = np.concatenate((cr_lines_mean, cr[-1:]), axis = 0)
new_cr[1::2,:] = cr_lines_mean
else:
cb_ratio = int(upsample_ratio[0]/upsample_ratio[1])
cr_ratio = int(upsample_ratio[0]/upsample_ratio[2])
cb_mean = (cb[:,:-1] + np.roll(cb, -1, 1)[:,:-1])/2
cb_mean = np.concatenate((cb_mean, cb[:,-1:]), axis = 1)
new_cb[:,::cb_ratio] = cb
new_cb[:,1::cb_ratio] = cb_mean
cr_mean = (cr[:,:-1] + np.roll(cr, -1, 1)[:,:-1])/2
cr_mean = np.concatenate((cr_mean, cr[:,-1:]), axis = 1)
new_cr[:,::cr_ratio] = cr
new_cr[:,1::cr_ratio] = cr_mean
return y,new_cb, new_cr
ratio = (4, 2, 2)
y_d, cb_d, cr_d = sub_sample(y, cb, cr, ratio)
show_image(y_d, cm_gray)
show_image(cb_d, cm_cb)
show_image(cr_d, cm_cr)
print("Y shape: " + str(y_d.shape),
"Cb shape: " + str(cb_d.shape),
"Cr shape: " + str(cr_d.shape),
sep = '\n'
)
Y shape: (304, 400) Cb shape: (304, 200) Cr shape: (304, 200)
ratio = (4, 2, 0)
y_d, cb_d, cr_d = sub_sample(y, cb, cr, ratio)
show_image(y_d, cm_gray)
show_image(cb_d, cm_cb)
show_image(cr_d, cm_cr)
print("Y shape: " + str(y_d.shape),
"Cb shape: " + str(cb_d.shape),
"Cr shape: " + str(cr_d.shape),
sep = '\n'
)
Y shape: (304, 400) Cb shape: (152, 200) Cr shape: (152, 200)
No downsampling 4:2:2, há uma redução do número de pixeis, nos canais Cb e Cr, de 2:1. No caso do 4:2:0, a taxa de compressão é maior, pois é de 4:1 nesses mesmos canais.
Esta etapa de downsampling é destrutiva, pelo que a contrapartida é que quanto mais alta for a compressão, mais alta é a taxa de destrutividade, pois quantos mais dados forem eliminados, menos fiel vai ser a reconstrução (tendo em conta os nosso preditores).
Assim, a escolha destas métricas deve ter em conta a qualidade que se pretende obter na reconstrução da imagem.
def dct_(a):
return dct(dct(a, norm="ortho").T, norm = "ortho").T
def idct_(a_dct):
return idct(idct(a_dct, norm='ortho').T, norm='ortho').T
# DCT aplicada na imagem toda
y_dct = dct_(y_d)
cb_dct = dct_(cb_d)
cr_dct = dct_(cr_d)
show_image(np.log(np.abs(y_dct) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dct) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dct) + 0.0001), cm_gray)
# Inverso da DCT
y_idct = idct_(y_dct)
cb_idct = idct_(cb_dct)
cr_idct = idct_(cr_dct)
show_image(y_idct, cm_gray)
show_image(cb_idct, cm_gray)
show_image(cr_idct, cm_gray)
# Reconstrução com copia da coluna da esquerda
y, cb, cr = const_up_sample(y_idct, cb_idct, cr_idct, ratio)
show_image(y, cm_gray)
show_image(cb, cm_gray)
show_image(cr, cm_gray)
#Recontrução com interpolação linear (média) da coluna anterior e seguinte
y, cb, cr = linear_up_sample(y_idct, cb_idct, cr_idct, ratio)
show_image(y, cm_gray)
show_image(cb, cm_gray)
show_image(cr, cm_gray)
A DCT, tal como a DFT, permite representar um sinal por valores de frequência com coeficientes associados. No entanto, a DCT tem a vantagem de conseguir concentrar muito mais a energia.
No caso, a DCT faz uma transformação do domínio do espaço para o domínio da frequência. Nesta transformação, os coeficientes das baixas frequências ficam presentes no canto superior esquerdo e os das altas frequências, no canto inferior direito.
Normalmente, no amplo espaço de uma imagem há vários tipos de transições, abruscas e suaves. Como é visivel em todas as imagens não há uma concentração de energia notória na DCT das imagens. Aliás, a matriz resultante demonstra-se semelhante ao ruído branco, que tem uma baixa taxa de compressão.
Semana 3¶
Aplicação da DCT em blocos 8x8 ou 64x64
- DCT em blocos e inverso
def dct_n_blocks(y, cb, cr, step):
dct_y = np.zeros(y.shape)
dct_cb = np.zeros(cb.shape)
dct_cr = np.zeros(cr.shape)
for i in range(0, y.shape[0], step):
for j in range(0, y.shape[1], step):
dct_y[i:i+step, j:j+step] = dct_(y[i:i+step, j:j+step])
for i in range(0, cb.shape[0], step):
for j in range(0, cb.shape[1], step):
dct_cb[i:i+step, j:j+step] = dct_(cb[i:i+step, j:j+step])
dct_cr[i:i+step, j:j+step] = dct_(cr[i:i+step, j:j+step])
return dct_y, dct_cb, dct_cr
def idct_n_blocks(y, cb, cr, step):
idct_y = np.zeros(y.shape)
idct_cb = np.zeros(cb.shape)
idct_cr = np.zeros(cr.shape)
for i in range(0, y.shape[0], step):
for j in range(0, y.shape[1], step):
idct_y[ i : i+step, j:j+step] = idct_(y[i:i+step, j:j+step])
for i in range(0, cb.shape[0], step):
for j in range(0, cb.shape[1], step):
idct_cb[ i : i+step, j : j+step] = idct_(cb[ i : i+step, j : j+step])
idct_cr[ i : i+step, j : j+step] = idct_(cr[ i : i+step, j : j+step])
return idct_y, idct_cb, idct_cr
# DCT em blocos
n = 64
y_dctn, cb_dctn, cr_dctn = dct_n_blocks(y, cb, cr, n)
show_image(np.log(np.abs(y_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn) + 0.0001), cm_gray)
# Inverter DCT em blocos
y_d, cb_d, cr_d = idct_n_blocks(y_dctn, cb_dctn, cr_dctn, n)
show_image(y_d, cm_gray)
show_image(cb_d, cm_gray)
show_image(cr_d, cm_gray)
# DCT em blocos
n = 8
y_dctn, cb_dctn, cr_dctn = dct_n_blocks(y, cb, cr, n)
show_image(np.log(np.abs(y_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn) + 0.0001), cm_gray)
# Inverter DCT em blocos
y_d, cb_d, cr_d = idct_n_blocks(y_dctn, cb_dctn, cr_dctn, n)
show_image(y_d, cm_gray)
show_image(cb_d, cm_gray)
show_image(cr_d, cm_gray)
A aplicação da DCT por bloco parte do princípio que num pequeno espaço da imagem, a variação vai ser muito baixa.
Ao contrário do que aconteceu no cenário anterior, em que foi aplicada a DCT em toda a imagem, no caso da aplicação da DCT em blocos 8x8 é clara a concentração da energia em cada bloco, o que aumenta o potencial de compressão.
No caso da transformada aplicada a blocos 64x64, também é visível a concentração da energia, mas não tanto como no caso 8x8. É visível, em cada bloco, algo semelhante a ruído branco.
Semana 4¶
Quantização da imagem e compressão dos coponentes DC da DCT com delta encoding
- Quantização dos canais Y, Cb e Cr com as matrizes adequadas
- delta encoding para as componentes DC da DCT.
- A compressão dos componentes AC não faz parte deste trabalho
def get_Qs():
Q_y = np.array(
[
[16, 11, 10, 16, 24, 40, 51, 61],
[12, 12, 14, 19, 26, 58, 60, 55],
[14, 13, 16, 24, 40, 57, 69, 56],
[14, 17, 22, 29, 51, 87, 80, 62],
[18, 22, 37, 56, 68, 109, 103, 77],
[24, 35, 55, 64, 81, 104, 113, 92],
[49, 64, 78, 87, 103, 121, 120, 101],
[72, 92, 95, 98, 112, 100, 103, 99],
]
)
Q_CbCr = np.array(
[
[17, 18, 24, 47, 99, 99, 99, 99],
[18, 21, 26, 66, 99, 99, 99, 99],
[24, 26, 56, 99, 99, 99, 99, 99],
[47, 66, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99],
[99, 99, 99, 99, 99, 99, 99, 99],
]
)
return Q_y, Q_CbCr
def quantitization(y, cb, cr, fator):
q_y, q_cbcr = get_Qs()
y_q = np.zeros(y.shape)
cb_q = np.zeros(cb.shape)
cr_q = np.zeros(cr.shape)
if fator>=50:
s = (100-fator)/50
else:
s = 50/fator
if s!=0:
q_y = np.round(s*q_y)
q_cbcr = np.round(s*q_cbcr)
else:
q_y[::,::] = 1
q_cbcr[::,::] = 1
q_y[q_y>255] = 255
q_cbcr[q_cbcr>255] = 255
for r in range(0,y.shape[0],8):
for c in range(0,y.shape[1],8):
y_q[r:r+8, c:c+8] = y[ r : r+8, c : c+8] / q_y
for r in range(0, cb.shape[0], 8):
for c in range(0,cb.shape[1], 8):
cb_q[ r : r+8, c : c+8] = cb[ r : r+8, c : c+8] / q_cbcr
cr_q[ r : r+8, c : c+8] = cr[ r : r+8, c : c+8] / q_cbcr
y_q = np.round(y_q)
cb_q = np.round(cb_q)
cr_q = np.round(cr_q)
return y_q, cb_q, cr_q
def quantitization_inverse(y, cb, cr, fator):
q_y, q_cbcr = get_Qs()
y_q = np.zeros(y.shape, dtype=float)
cb_q = np.zeros(cb.shape, dtype= float)
cr_q = np.zeros(cr.shape, dtype=float)
if fator>=50:
s = (100-fator)/50
else:
s = 50/fator
if s!=0:
q_y = np.round(s*q_y)
q_cbcr = np.round(s*q_cbcr)
else:
q_y[::,::] = 1
q_cbcr[::,::] = 1
q_y[ q_y > 255 ] = 255
q_cbcr[ q_cbcr > 255 ] = 255
for r in range(0, y.shape[0], 8):
for c in range(0, y.shape[1], 8):
y_q[ r : r+8, c : c+8] = y[ r : r+8, c : c+8]*q_y
for r in range(0, cb.shape[0], 8):
for c in range(0, cb.shape[1], 8):
cb_q[ r : r+8, c : c+8] = cb[ r : r+8, c : c+8] * q_cbcr
cr_q[ r : r+8, c : c+8] = cr[ r : r+8, c : c+8] * q_cbcr
y_q = np.round(y_q)
cb_q = np.round(cb_q)
cr_q = np.round(cr_q)
return y_q, cb_q, cr_q
# Fator de qualidade para a quantização da imagem
fator = 100
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)
show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
# Fator de qualidade para a quantização da imagem
fator = 75
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)
show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
# Fator de qualidade para a quantização da imagem
fator = 50
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)
show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
# Fator de qualidade para a quantização da imagem
fator = 25
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)
show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
# Fator de qualidade para a quantização da imagem
fator = 10
y_dctn_q, cb_dctn_q, cr_dctn_q = quantitization(y_dctn, cb_dctn, cr_dctn, fator)
y_dctn_qi, cb_dctn_qi, cr_dctn_qi = quantitization_inverse(y_dctn_q, cb_dctn_q, cr_dctn_q, fator)
show_image(np.log(np.abs(y_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cb_dctn_q) + 0.0001), cm_gray)
show_image(np.log(np.abs(cr_dctn_q) + 0.0001), cm_gray)
É possível verificar que à medida que o fator de qualidade diminuiu, também diminuiu o número de frequências que reprentam cada bloco. Em casos em que o fator é reduzido há um grande número de blocos que aparentam ser representados apenas pelo seu componente DC.
Com a imagem representada por um menor número de diferentes valores, o potencial de compressão aumenta, pois a distribuição dos valores diminuiu, tal como a entropia. Assim, quanto menor for o fator de qualidade, maior é o potencial de compressão.
A crominância tem menos relevância na perceção da imagem ao olho humano, por isso os canais Cb e Cr podem sofrer mais perda de informação, como indicado anteriormente.
É visível que há proveito desta característica humana uma vez que em todos os exemplos o canal Y mantém mais informação que os canais Cb e Cr. Há uma tentativa de minimizar a perda de informação relevante sobre a luminância na compressão destrutiva.
# Delta encoding aplicado aos coeficientes DC da DCT em blocos
def diff(bloco):
x, y = bloco.shape
bloco = bloco.flatten().astype(float)
bloco[1:] -= bloco[:-1]
return bloco.reshape(x,y)
def DPCM(y,cb,cr):
yBlocks = y.copy()
cbBlocks = cb.copy()
crBlocks = cr.copy()
yBlocks[::8,::8] = diff(y[::8,::8])
cbBlocks[::8,::8] = diff(cb[::8,::8])
crBlocks[::8,::8] = diff(cr[::8,::8])
return yBlocks, cbBlocks, crBlocks
def diff_reverse(bloco):
x, y = bloco.shape
bloco = bloco.flatten().astype(float)
for i in range(1,bloco.shape[0]):
bloco[i] = bloco[i] + bloco[i-1]
return bloco.reshape(x,y)
def DPCM_reverse(yBlocks, cbBlocks, crBlocks):
y = yBlocks.copy()
cr = crBlocks.copy()
cb = cbBlocks.copy()
y[::8,::8] = diff_reverse(yBlocks[::8,::8])
cr[::8,::8] = diff_reverse(crBlocks[::8,::8])
cb[::8,::8] = diff_reverse(cbBlocks[::8,::8])
return y, cb, cr
d_y, d_cb, d_cr = DPCM(y_dctn_q, cb_dctn_q, cr_dctn_q)
id_y, id_cb, id_cr = DPCM_reverse(d_y, d_cb, d_cr)
show_image(np.log(np.abs(d_y) + 0.0001), cm_gray)
show_image(np.log(np.abs(d_cb) + 0.0001), cm_gray)
show_image(np.log(np.abs(d_cr) + 0.0001), cm_gray)
Após a aplicação da codificação DPCM, o número de blocos contíguos totalmente nulos aumentou bastante. Isto potencializa a compressão de algorítmos como o RLE.
De notar que esta compressão não é destrutiva, pelo que este passo ajuda a compressão sem perder qualidade na reconstrução da imagem, através de um método simples de aplicar.
Semana 5¶
Aplicação de todas as funções criadas antriormente com os fatores de qualidade 10, 25, 50, 75, 100
- Cálculo das métricas de distorção (MSE, RMSE, SNR, PSNR) para todas as imagens e fatores de qualidade
def MSE(orig, compressed):
return np.sum(np.square(orig-compressed))/np.prod(orig.shape) * 3
def RMSE(orig, compressed):
return np.sqrt(MSE(orig,compressed))
def SNR(orig, compressed):
P = np.sum(np.square(orig))/np.prod(orig.shape)
return 10*np.log10(P/MSE(orig, compressed))
def PSNR(orig, compressed):
return 10*np.log10(np.square(orig).max()/MSE(orig,compressed))
def encoder(filename, ratio, n, factor):
img = read_image(filename)
shape, pad_img = padding(img)
ycbcr = rgb_to_ycbcr(pad_img)
y, cb, cr = separate_components(ycbcr)
y_d, cb_d, cr_d = sub_sample(y, cb, cr, ratio)
y_dct, cb_dct, cr_dct = dct_n_blocks(y_d, cb_d, cr_d, n)
y_dct_q, cb_dct_q, cr_dct_q = quantitization(y_dct, cb_dct, cr_dct, factor)
delta_y, delta_cb, delta_cr = DPCM(y_dct_q, cb_dct_q, cr_dct_q)
return delta_y, delta_cb, delta_cr, shape
def decoder(delta_y, delta_cb, delta_cr, ratio, shape, n, factor):
y_dct_q, cb_dct_q, cr_dct_q = DPCM_reverse(delta_y, delta_cb, delta_cr)
y_dct, cb_dct, cr_dct = quantitization_inverse(y_dct_q, cb_dct_q, cr_dct_q, factor)
y_idct, cb_idct, cr_idct = idct_n_blocks(y_dct, cb_dct, cr_dct, n) # DCT em blocos n x n
#y, cb, cr = const_up_sample(y_idct, cb_idct, cr_idct, ratio) # Upsample: cópia da coluna á esquerda
y, cb, cr = linear_up_sample(y_idct, cb_idct, cr_idct, ratio) # Upsample: média entre o anterior e o seguinte
ycbcr = join_components(y, cb, cr)
img_rgb = ycbcr_to_rgb(ycbcr)
img = unpad(img_rgb, shape)
return img
%matplotlib inline
f = "imagens/barn_mountains.bmp"
ratio = (4,2,2)
n = 8
fator = 75
x,y,z,shape = encoder(f,ratio,n,fator)
img = decoder(x,y,z,ratio,shape,n,fator)
show_image(img)
print("MSE:", MSE(read_image(f).astype(float), img.astype(float)),
"RMSE:", RMSE(read_image(f).astype(float), img.astype(float)),
"SNR:", SNR(read_image(f).astype(float), img.astype(float)),
"PSNR:",PSNR(read_image(f).astype(float), img.astype(float)),
sep = "\n")
MSE: 152.02205387205387 RMSE: 12.329722376114308 SNR: 20.595997303892048 PSNR: 26.311737651588494
imagens = ['imagens/barn_mountains.bmp', 'imagens/peppers.bmp', 'imagens/logo.bmp']
ratio = (4,2,2)
n = 8
fatores = [10, 25, 50, 75, 100]
for file in imagens:
print(f"{file.split('/')[1]}")
f, ax = plt.subplots(nrows = 1, ncols = 5, figsize = (20, 5))
f.suptitle(f"{file.split('/')[1]}", fontsize = 16)
df = pd.DataFrame(index = fatores, columns = ["MSE", "RMSE", "SNR", "PSNR"])
for index, fator in enumerate(fatores):
y, cb, cr, shape = encoder(file, ratio, n, fator)
i_jpeg = decoder(y, cb, cr, ratio, shape, n, fator)
ycbcr = rgb_to_ycbcr(i_jpeg)
y, cb, cr = separate_components(ycbcr)
img = read_image(file)
ycbcr = rgb_to_ycbcr(img)
y_orig, cb_orig, cr_orig = separate_components(ycbcr)
y_diff = np.abs(y_orig - y)
ax[index].title.set_text(f"{file.split('/')[1]}: {fator}")
y_diff[-1,-1] = 255
y_diff[-1, -2] = 0
ax[index].imshow(y_diff, cm_gray)
df["MSE"][fator] = MSE(img.astype(float), i_jpeg.astype(float))
df["RMSE"][fator] = RMSE(img.astype(float), i_jpeg.astype(float))
df["SNR"][fator] = SNR(img.astype(float), i_jpeg.astype(float))
df["PSNR"][fator] = PSNR(img.astype(float), i_jpeg.astype(float))
print(df)
for i in ax:
i.set_axis_off()
barn_mountains.bmp
MSE RMSE SNR PSNR
10 707.49287 26.598738 13.917843 19.633583
25 398.153131 19.953775 16.414562 22.130302
50 261.369487 16.166926 18.242514 23.958255
75 152.022054 12.329722 20.595997 26.311738
100 11.170572 3.342241 31.934309 37.650049
peppers.bmp
MSE RMSE SNR PSNR
10 283.115982 16.826051 15.635096 23.61116
25 127.319875 11.283611 19.105778 27.081842
50 77.324849 8.793455 21.271549 29.247613
75 50.02124 7.07257 23.163195 31.139259
100 6.346832 2.519292 32.129169 40.105233
logo.bmp
MSE RMSE SNR PSNR
10 157.304726 12.542118 24.593488 26.163386
25 65.218448 8.075794 28.417201 29.987099
50 43.806448 6.618644 30.145526 31.715423
75 23.667495 4.864925 32.819383 34.389281
100 6.184128 2.486791 38.648121 40.218019
Tanto o ruído quanto a diferença entre os pixeis da imagem original e da imagem reconstruida serão tanto maiores quanto menor for fator de qualidade. Desta forma, pode-se concluir que o SNR e o PSNR terão valores superiores em imagens cujo fator de qualidade seja alto e o MSE e RMSE terão valores tanto maiores quanto menor for o fator de qualidade.